//
// Saru.Base.Math.matrix3.cs: 3 Dimensional matrices and operations
//
// Author:
//   Jan Van Sweevelt (vansweej@gmail.com)
//
// Copyright (C) 2007
// 

using System;

namespace Saru.Base.Math {

    public class Matrix4f {

	// Constructors

        // create identity matrix
	public Matrix4f() {
	    tuple = new float[16];
	    MakeIdentity();
	}

	// copy constructor
	public Matrix4f(Matrix4f M) {
	    tuple = new float[16];
	    tuple[0] = M.tuple[0];
	    tuple[1] = M.tuple[1];
	    tuple[2] = M.tuple[2];
	    tuple[3] = M.tuple[3];
	    tuple[4] = M.tuple[4];
	    tuple[5] = M.tuple[5];
	    tuple[6] = M.tuple[6];
	    tuple[7] = M.tuple[7];
	    tuple[8] = M.tuple[8];
	    tuple[9] = M.tuple[9];
	    tuple[10] = M.tuple[10];
	    tuple[11] = M.tuple[11];
	    tuple[12] = M.tuple[12];
	    tuple[13] = M.tuple[13];
	    tuple[14] = M.tuple[14];
	    tuple[15] = M.tuple[15];
	}

	// construct from floating points
	public Matrix4f(float M00, float M01, float M02, float M03, float M10, float M11, float M12, float M13, float M20, float M21, float M22, float M23, float M30, float M31, float M32, float M33) {
	    tuple = new float[16];
	    tuple[0] = M00;
	    tuple[1] = M01;
	    tuple[2] = M02;
	    tuple[3] = M03;
	    tuple[4] = M10;
	    tuple[5] = M11;
	    tuple[6] = M12;
	    tuple[7] = M13;
	    tuple[8] = M20;
	    tuple[9] = M21;
	    tuple[10] = M22;
	    tuple[11] = M23;
	    tuple[12] = M30;
	    tuple[13] = M31;
	    tuple[14] = M32;
	    tuple[15] = M33;
	}

	// construct from vectors
	// colums means:
	// true: vectors are colums of the matrix
	// false: vectors are rows of the matrix
	public Matrix4f(Vector4f V1, Vector4f V2, Vector4f V3, Vector4f V4, bool colums) {
	    tuple = new float[16];
	    if(colums) {
		tuple[0] = V1.X;
		tuple[1] = V2.X;
		tuple[2] = V3.X;
		tuple[3] = V4.X;
		tuple[4] = V1.Y;
		tuple[5] = V2.Y;
		tuple[6] = V3.Y;
		tuple[7] = V4.Y;
		tuple[8] = V1.Z;		
		tuple[9] = V2.Z;
		tuple[10] = V3.Z;
		tuple[11] = V4.Z;
		tuple[12] = V1.W;
		tuple[13] = V2.W;
		tuple[14] = V3.W;
		tuple[15] = V4.W;
	    }
	    else {
		tuple[0] = V1.X;
		tuple[1] = V1.Y;
		tuple[2] = V1.Z;
		tuple[3] = V1.W;		
		tuple[4] = V2.X;
		tuple[5] = V2.Y;
		tuple[6] = V2.Z;
		tuple[7] = V2.W;
		tuple[8] = V3.X;		
		tuple[9] = V3.Y;
		tuple[10] = V3.Z;
		tuple[11] = V3.W;
		tuple[12] = V4.X;		
		tuple[13] = V4.Y;
		tuple[14] = V4.Z;
		tuple[15] = V4.W;
	    }
	}

	// create a diagonal matrix
	public Matrix4f(float M00, float M11, float M22, float M33) {
	    tuple = new float[16];

	    MakeDiagonal(M00, M11, M22, M33);
	}

	public void MakeZero() {
	    tuple[0] = 0.0f;
	    tuple[1] = 0.0f;
	    tuple[2] = 0.0f;
	    tuple[3] = 0.0f;
	    tuple[4] = 0.0f;
	    tuple[5] = 0.0f;
	    tuple[6] = 0.0f;
	    tuple[7] = 0.0f;
	    tuple[8] = 0.0f;
	    tuple[9] = 0.0f;
	    tuple[10] = 0.0f;
	    tuple[11] = 0.0f;
	    tuple[12] = 0.0f;
	    tuple[13] = 0.0f;
	    tuple[14] = 0.0f;
	    tuple[15] = 0.0f;
	}

	public void MakeIdentity() {
	    tuple[0] = 1.0f;
	    tuple[1] = 0.0f;
	    tuple[2] = 0.0f;
	    tuple[3] = 0.0f;
	    tuple[4] = 0.0f;
	    tuple[5] = 1.0f;
	    tuple[6] = 0.0f;
	    tuple[7] = 0.0f;
	    tuple[8] = 0.0f;
	    tuple[9] = 0.0f;
	    tuple[10] = 1.0f;
	    tuple[11] = 0.0f;
	    tuple[12] = 0.0f;
	    tuple[13] = 0.0f;
	    tuple[14] = 0.0f;
	    tuple[15] = 1.0f;
	}

	public void MakeDiagonal(float M00, float M11, float M22, float M33) {
	    tuple[0] = M00;
	    tuple[1] = 0.0f;
	    tuple[2] = 0.0f;
	    tuple[3] = 0.0f;
	    tuple[4] = 0.0f;
	    tuple[5] = M11;
	    tuple[6] = 0.0f;
	    tuple[7] = 0.0f;
	    tuple[8] = 0.0f;
	    tuple[9] = 0.0f;
	    tuple[10] = M22;
	    tuple[11] = 0.0f;
	    tuple[12] = 0.0f;
	    tuple[13] = 0.0f;
	    tuple[14] = 0.0f;
	    tuple[15] = M33;

	}

       // accessors
	public float[] Matrix4 {
	    get {
		return tuple;
	    }

	    set {
		tuple = value;
	    }
	}

	public float Get(short Column, short Row) {
	    return tuple[Column + 4 * Row];
	}

	public void Set(float Value, short Column, short Row) {
	    tuple[Column + 4 * Row] = Value;
	}

	public Vector4f GetRow(short Row) {
	    int i0 = 4*Row;
	    return new Vector4f(tuple[i0], tuple[i0+1], tuple[i0+2], tuple[i0+3]);
	}

	public void SetRow(short Row, Vector4f V) {
	    int i0 = 4*Row;
	    tuple[i0] = V.X;
	    tuple[i0+1] = V.Y;
	    tuple[i0+2] = V.Z;
	    tuple[i0+3] = V.W;
	}

	// comparison
	public static bool operator==(Matrix4f M1, Matrix4f M2) {
	    return M1.tuple == M2.tuple;
	}

	public static bool operator!=(Matrix4f M1, Matrix4f M2) {
	    return M1.tuple != M2.tuple;
	}

	public static bool operator<(Matrix4f M1, Matrix4f M2) {
	    return (M1.tuple[0] < M2.tuple[0] &&
		    M1.tuple[1] < M2.tuple[1] &&
		    M1.tuple[2] < M2.tuple[2] &&
		    M1.tuple[3] < M2.tuple[3] &&
		    M1.tuple[4] < M2.tuple[4] &&
		    M1.tuple[5] < M2.tuple[5] &&
		    M1.tuple[6] < M2.tuple[6] &&
		    M1.tuple[7] < M2.tuple[7] &&
		    M1.tuple[8] < M2.tuple[8] &&
		    M1.tuple[9] < M2.tuple[9] &&
		    M1.tuple[10] < M2.tuple[10] &&
		    M1.tuple[11] < M2.tuple[11] &&
		    M1.tuple[12] < M2.tuple[12] &&
		    M1.tuple[13] < M2.tuple[13] &&
		    M1.tuple[14] < M2.tuple[14] &&
		    M1.tuple[15] < M2.tuple[15]);
	}

	public static bool operator>(Matrix4f M1, Matrix4f M2) {
	    return (M1.tuple[0] > M2.tuple[0] &&
		    M1.tuple[1] > M2.tuple[1] &&
		    M1.tuple[2] > M2.tuple[2] &&
		    M1.tuple[3] > M2.tuple[3] &&
		    M1.tuple[4] > M2.tuple[4] &&
		    M1.tuple[5] > M2.tuple[5] &&
		    M1.tuple[6] > M2.tuple[6] &&
		    M1.tuple[7] > M2.tuple[7] &&
		    M1.tuple[8] > M2.tuple[8] &&
		    M1.tuple[9] > M2.tuple[9] &&
		    M1.tuple[10] > M2.tuple[10] &&
		    M1.tuple[11] > M2.tuple[11] &&
		    M1.tuple[12] > M2.tuple[12] &&
		    M1.tuple[13] > M2.tuple[13] &&
		    M1.tuple[14] > M2.tuple[14] &&
		    M1.tuple[15] > M2.tuple[15]);
	}

	public static bool operator<=(Matrix4f M1, Matrix4f M2) {
	    return (M1.tuple[0] <= M2.tuple[0] &&
		    M1.tuple[1] <= M2.tuple[1] &&
		    M1.tuple[2] <= M2.tuple[2] &&
		    M1.tuple[3] <= M2.tuple[3] &&
		    M1.tuple[4] <= M2.tuple[4] &&
		    M1.tuple[5] <= M2.tuple[5] &&
		    M1.tuple[6] <= M2.tuple[6] &&
		    M1.tuple[7] <= M2.tuple[7] &&
		    M1.tuple[8] <= M2.tuple[8] &&
		    M1.tuple[9] <= M2.tuple[9] &&
		    M1.tuple[10] <= M2.tuple[10] &&
		    M1.tuple[11] <= M2.tuple[11] &&
		    M1.tuple[12] <= M2.tuple[12] &&
		    M1.tuple[13] <= M2.tuple[13] &&
		    M1.tuple[14] <= M2.tuple[14] &&
		    M1.tuple[15] <= M2.tuple[15]);
	}

	public static bool operator>=(Matrix4f M1, Matrix4f M2) {
	    return (M1.tuple[0] >= M2.tuple[0] &&
		    M1.tuple[1] >= M2.tuple[1] &&
		    M1.tuple[2] >= M2.tuple[2] &&
		    M1.tuple[3] >= M2.tuple[3] &&
		    M1.tuple[4] >= M2.tuple[4] &&
		    M1.tuple[5] >= M2.tuple[5] &&
		    M1.tuple[6] >= M2.tuple[6] &&
		    M1.tuple[7] >= M2.tuple[7] &&
		    M1.tuple[8] >= M2.tuple[8] &&
		    M1.tuple[9] >= M2.tuple[9] &&
		    M1.tuple[10] >= M2.tuple[10] &&
		    M1.tuple[11] >= M2.tuple[11] &&
		    M1.tuple[12] >= M2.tuple[12] &&
		    M1.tuple[13] >= M2.tuple[13] &&
		    M1.tuple[14] >= M2.tuple[14] &&
		    M1.tuple[15] >= M2.tuple[15]);
	}

	public override bool Equals(object value) {
	    if(value == null) return false;
	    if (GetType() == value.GetType()) return this == (Matrix4f)value;
	    
	    return false;
	}

	public override int GetHashCode() {
	    return tuple.GetHashCode();
	}

	// operators
	public static Matrix4f operator+(Matrix4f M1, Matrix4f M2) {
	    return new Matrix4f(M1.tuple[0] + M2.tuple[0],
				M1.tuple[1] + M2.tuple[1],
				M1.tuple[2] + M2.tuple[2],
				M1.tuple[3] + M2.tuple[3],
				M1.tuple[4] + M2.tuple[4],
				M1.tuple[5] + M2.tuple[5],
				M1.tuple[6] + M2.tuple[6],
				M1.tuple[7] + M2.tuple[7],
				M1.tuple[8] + M2.tuple[8],
				M1.tuple[9] + M2.tuple[9],
				M1.tuple[10] + M2.tuple[10],
				M1.tuple[11] + M2.tuple[11],
				M1.tuple[12] + M2.tuple[12],
				M1.tuple[13] + M2.tuple[13],
				M1.tuple[14] + M2.tuple[14],
				M1.tuple[15] + M2.tuple[15]);
	}

	public static Matrix4f operator-(Matrix4f M1, Matrix4f M2) {
	    return new Matrix4f(M1.tuple[0] - M2.tuple[0],
				M1.tuple[1] - M2.tuple[1],
				M1.tuple[2] - M2.tuple[2],
				M1.tuple[3] - M2.tuple[3],
				M1.tuple[4] - M2.tuple[4],
				M1.tuple[5] - M2.tuple[5],
				M1.tuple[6] - M2.tuple[6],
				M1.tuple[7] - M2.tuple[7],
				M1.tuple[8] - M2.tuple[8],
				M1.tuple[9] - M2.tuple[9],
				M1.tuple[10] - M2.tuple[10],
				M1.tuple[11] - M2.tuple[11],
				M1.tuple[12] - M2.tuple[12],
				M1.tuple[13] - M2.tuple[13],
				M1.tuple[14] - M2.tuple[14],
				M1.tuple[15] - M2.tuple[15]);
	}

	public static Matrix4f operator*(Matrix4f M1, Matrix4f M2) {
	    return new Matrix4f(M1.tuple[0] * M2.tuple[0] + 
				M1.tuple[1] * M2.tuple[4] + 
				M1.tuple[2] * M2.tuple[8] +
				M1.tuple[3] * M2.tuple[12],

				M1.tuple[0] * M2.tuple[1] + 
				M1.tuple[1] * M2.tuple[5] + 
				M1.tuple[2] * M2.tuple[9] +
				M1.tuple[3] * M2.tuple[13],

				M1.tuple[0] * M2.tuple[2] + 
				M1.tuple[1] * M2.tuple[6] + 
				M1.tuple[2] * M2.tuple[10] +
				M1.tuple[3] * M2.tuple[14],

				M1.tuple[0] * M2.tuple[3] + 
				M1.tuple[1] * M2.tuple[7] + 
				M1.tuple[2] * M2.tuple[11] +
				M1.tuple[3] * M2.tuple[15],

				M1.tuple[4] * M2.tuple[0] + 
				M1.tuple[5] * M2.tuple[4] + 
				M1.tuple[6] * M2.tuple[8] +
				M1.tuple[7] * M2.tuple[12],

				M1.tuple[4] * M2.tuple[1] + 
				M1.tuple[5] * M2.tuple[5] + 
				M1.tuple[6] * M2.tuple[9] +
				M1.tuple[7] * M2.tuple[13],

				M1.tuple[4] * M2.tuple[2] + 
				M1.tuple[5] * M2.tuple[6] + 
				M1.tuple[6] * M2.tuple[10] +
				M1.tuple[7] * M2.tuple[14],

				M1.tuple[4] * M2.tuple[3] + 
				M1.tuple[5] * M2.tuple[7] + 
				M1.tuple[6] * M2.tuple[11] +
				M1.tuple[7] * M2.tuple[15],

				M1.tuple[8] * M2.tuple[0] + 
				M1.tuple[9] * M2.tuple[4] + 
				M1.tuple[10] * M2.tuple[8] +
				M1.tuple[11] * M2.tuple[12],

				M1.tuple[8] * M2.tuple[1] + 
				M1.tuple[9] * M2.tuple[5] + 
				M1.tuple[10] * M2.tuple[9] +
				M1.tuple[11] * M2.tuple[13],

				M1.tuple[8] * M2.tuple[2] + 
				M1.tuple[9] * M2.tuple[6] + 
				M1.tuple[10] * M2.tuple[10] +
				M1.tuple[11] * M2.tuple[14],

				M1.tuple[8] * M2.tuple[3] + 
				M1.tuple[9] * M2.tuple[7] + 
				M1.tuple[10] * M2.tuple[11] +
				M1.tuple[11] * M2.tuple[15],

				M1.tuple[12] * M2.tuple[0] + 
				M1.tuple[13] * M2.tuple[4] + 
				M1.tuple[14] * M2.tuple[8] +
				M1.tuple[15] * M2.tuple[12],

				M1.tuple[12] * M2.tuple[1] + 
				M1.tuple[13] * M2.tuple[5] + 
				M1.tuple[14] * M2.tuple[9] +
				M1.tuple[15] * M2.tuple[13],

				M1.tuple[12] * M2.tuple[2] + 
				M1.tuple[13] * M2.tuple[6] + 
				M1.tuple[14] * M2.tuple[10] +
				M1.tuple[15] * M2.tuple[14],

				M1.tuple[12] * M2.tuple[3] + 
				M1.tuple[13] * M2.tuple[7] + 
				M1.tuple[14] * M2.tuple[11] +
				M1.tuple[15] * M2.tuple[15]);
	}

	public static Matrix4f operator*(Matrix4f M, float Scalar) {
	    return new Matrix4f(M.tuple[0] * Scalar,
				M.tuple[1] * Scalar,
				M.tuple[2] * Scalar,
				M.tuple[3] * Scalar,
				M.tuple[4] * Scalar,
				M.tuple[5] * Scalar,
				M.tuple[6] * Scalar,
				M.tuple[7] * Scalar,
				M.tuple[8] * Scalar,
				M.tuple[9] * Scalar,
				M.tuple[10] * Scalar,
				M.tuple[11] * Scalar,
				M.tuple[12] * Scalar,
				M.tuple[13] * Scalar,
				M.tuple[14] * Scalar,
				M.tuple[15] * Scalar);
	}

	public static Matrix4f operator/(Matrix4f M, float Scalar) {
	    if(Scalar != 0.0f) {
		float reciproke = 1.0f/Scalar;

		return new Matrix4f(M.tuple[0] * reciproke,
				    M.tuple[1] * reciproke,
				    M.tuple[2] * reciproke,
				    M.tuple[3] * reciproke,
				    M.tuple[4] * reciproke,
				    M.tuple[5] * reciproke,
				    M.tuple[6] * reciproke,
				    M.tuple[7] * reciproke,
				    M.tuple[8] * reciproke,
				    M.tuple[9] * reciproke,
				    M.tuple[10] * reciproke,
				    M.tuple[11] * reciproke,
				    M.tuple[12] * reciproke,
				    M.tuple[13] * reciproke,
				    M.tuple[14] * reciproke,
				    M.tuple[15] * reciproke);
	    }
	    else {
		return new Matrix4f(float.MaxValue,
				    float.MaxValue,
				    float.MaxValue,
				    float.MaxValue,
				    float.MaxValue,
				    float.MaxValue,
				    float.MaxValue,
				    float.MaxValue,
				    float.MaxValue,
				    float.MaxValue,
				    float.MaxValue,
				    float.MaxValue,
				    float.MaxValue,
				    float.MaxValue,
				    float.MaxValue,
				    float.MaxValue);
	    }
	}

	public static Vector4f operator*(Matrix4f M, Vector4f V) {
	    return new Vector4f(M.tuple[0] * V.X + 
				M.tuple[1] * V.Y + 
				M.tuple[2] * V.Z +
				M.tuple[3] * V.W,
				
				M.tuple[4] * V.X + 
				M.tuple[5] * V.Y + 
				M.tuple[6] * V.Z +
				M.tuple[7] * V.W,
				
				M.tuple[8] * V.X + 
				M.tuple[9] * V.Y + 
				M.tuple[10] * V.Z +
				M.tuple[11] * V.W,

				M.tuple[12] * V.X + 
				M.tuple[13] * V.Y + 
				M.tuple[14] * V.Z +
				M.tuple[15] * V.W);
	}

	public static Matrix4f operator-(Matrix4f M) {
	    return new Matrix4f(-M.tuple[0],
				-M.tuple[1],
				-M.tuple[2],
				-M.tuple[3],
				-M.tuple[4],
				-M.tuple[5],
				-M.tuple[6],
				-M.tuple[7],
				-M.tuple[8],
				-M.tuple[9],
				-M.tuple[10],
				-M.tuple[11],
				-M.tuple[12],
				-M.tuple[13],
				-M.tuple[14],
				-M.tuple[15]);
	}

	public Matrix4f Transpose() {
	    return new Matrix4f(tuple[0],
				tuple[4],
				tuple[8],
				tuple[12],
				tuple[1],
				tuple[5],
				tuple[9],
				tuple[13],
				tuple[2],
				tuple[6],
				tuple[10],
				tuple[14],
				tuple[3],
				tuple[7],
				tuple[11],
				tuple[15]);
	}

	private float[] tuple;
    }

    public class Matrix4d {

	// Constructors

        // create identity matrix
	public Matrix4d() {
	    tuple = new double[16];
	    MakeIdentity();
	}

	// copy constructor
	public Matrix4d(Matrix4d M) {
	    tuple = new double[16];
	    tuple[0] = M.tuple[0];
	    tuple[1] = M.tuple[1];
	    tuple[2] = M.tuple[2];
	    tuple[3] = M.tuple[3];
	    tuple[4] = M.tuple[4];
	    tuple[5] = M.tuple[5];
	    tuple[6] = M.tuple[6];
	    tuple[7] = M.tuple[7];
	    tuple[8] = M.tuple[8];
	    tuple[9] = M.tuple[9];
	    tuple[10] = M.tuple[10];
	    tuple[11] = M.tuple[11];
	    tuple[12] = M.tuple[12];
	    tuple[13] = M.tuple[13];
	    tuple[14] = M.tuple[14];
	    tuple[15] = M.tuple[15];
	}

	// construct from doubleing points
	public Matrix4d(double M00, double M01, double M02, double M03, double M10, double M11, double M12, double M13, double M20, double M21, double M22, double M23, double M30, double M31, double M32, double M33) {
	    tuple = new double[16];
	    tuple[0] = M00;
	    tuple[1] = M01;
	    tuple[2] = M02;
	    tuple[3] = M03;
	    tuple[4] = M10;
	    tuple[5] = M11;
	    tuple[6] = M12;
	    tuple[7] = M13;
	    tuple[8] = M20;
	    tuple[9] = M21;
	    tuple[10] = M22;
	    tuple[11] = M23;
	    tuple[12] = M30;
	    tuple[13] = M31;
	    tuple[14] = M32;
	    tuple[15] = M33;
	}

	// construct from vectors
	// colums means:
	// true: vectors are colums of the matrix
	// false: vectors are rows of the matrix
	public Matrix4d(Vector4d V1, Vector4d V2, Vector4d V3, Vector4d V4, bool colums) {
	    tuple = new double[16];
	    if(colums) {
		tuple[0] = V1.X;
		tuple[1] = V2.X;
		tuple[2] = V3.X;
		tuple[3] = V4.X;
		tuple[4] = V1.Y;
		tuple[5] = V2.Y;
		tuple[6] = V3.Y;
		tuple[7] = V4.Y;
		tuple[8] = V1.Z;		
		tuple[9] = V2.Z;
		tuple[10] = V3.Z;
		tuple[11] = V4.Z;
		tuple[12] = V1.W;
		tuple[13] = V2.W;
		tuple[14] = V3.W;
		tuple[15] = V4.W;
	    }
	    else {
		tuple[0] = V1.X;
		tuple[1] = V1.Y;
		tuple[2] = V1.Z;
		tuple[3] = V1.W;		
		tuple[4] = V2.X;
		tuple[5] = V2.Y;
		tuple[6] = V2.Z;
		tuple[7] = V2.W;
		tuple[8] = V3.X;		
		tuple[9] = V3.Y;
		tuple[10] = V3.Z;
		tuple[11] = V3.W;
		tuple[12] = V4.X;		
		tuple[13] = V4.Y;
		tuple[14] = V4.Z;
		tuple[15] = V4.W;
	    }
	}

	// create a diagonal matrix
	public Matrix4d(double M00, double M11, double M22, double M33) {
	    tuple = new double[16];

	    MakeDiagonal(M00, M11, M22, M33);
	}

	public void MakeZero() {
	    tuple[0] = 0.0;
	    tuple[1] = 0.0;
	    tuple[2] = 0.0;
	    tuple[3] = 0.0;
	    tuple[4] = 0.0;
	    tuple[5] = 0.0;
	    tuple[6] = 0.0;
	    tuple[7] = 0.0;
	    tuple[8] = 0.0;
	    tuple[9] = 0.0;
	    tuple[10] = 0.0;
	    tuple[11] = 0.0;
	    tuple[12] = 0.0;
	    tuple[13] = 0.0;
	    tuple[14] = 0.0;
	    tuple[15] = 0.0;
	}

	public void MakeIdentity() {
	    tuple[0] = 1.0;
	    tuple[1] = 0.0;
	    tuple[2] = 0.0;
	    tuple[3] = 0.0;
	    tuple[4] = 0.0;
	    tuple[5] = 1.0;
	    tuple[6] = 0.0;
	    tuple[7] = 0.0;
	    tuple[8] = 0.0;
	    tuple[9] = 0.0;
	    tuple[10] = 1.0;
	    tuple[11] = 0.0;
	    tuple[12] = 0.0;
	    tuple[13] = 0.0;
	    tuple[14] = 0.0;
	    tuple[15] = 1.0;
	}

	public void MakeDiagonal(double M00, double M11, double M22, double M33) {
	    tuple[0] = M00;
	    tuple[1] = 0.0;
	    tuple[2] = 0.0;
	    tuple[3] = 0.0;
	    tuple[4] = 0.0;
	    tuple[5] = M11;
	    tuple[6] = 0.0;
	    tuple[7] = 0.0;
	    tuple[8] = 0.0;
	    tuple[9] = 0.0;
	    tuple[10] = M22;
	    tuple[11] = 0.0;
	    tuple[12] = 0.0;
	    tuple[13] = 0.0;
	    tuple[14] = 0.0;
	    tuple[15] = M33;

	}

       // accessors
	public double[] Matrix4 {
	    get {
		return tuple;
	    }

	    set {
		tuple = value;
	    }
	}

	public double Get(short Column, short Row) {
	    return tuple[Column + 4 * Row];
	}

	public void Set(double Value, short Column, short Row) {
	    tuple[Column + 4 * Row] = Value;
	}

	public Vector4d GetRow(short Row) {
	    int i0 = 4*Row;
	    return new Vector4d(tuple[i0], tuple[i0+1], tuple[i0+2], tuple[i0+3]);
	}

	public void SetRow(short Row, Vector4d V) {
	    int i0 = 4*Row;
	    tuple[i0] = V.X;
	    tuple[i0+1] = V.Y;
	    tuple[i0+2] = V.Z;
	    tuple[i0+3] = V.W;
	}

	// comparison
	public static bool operator==(Matrix4d M1, Matrix4d M2) {
	    return M1.tuple == M2.tuple;
	}

	public static bool operator!=(Matrix4d M1, Matrix4d M2) {
	    return M1.tuple != M2.tuple;
	}

	public static bool operator<(Matrix4d M1, Matrix4d M2) {
	    return (M1.tuple[0] < M2.tuple[0] &&
		    M1.tuple[1] < M2.tuple[1] &&
		    M1.tuple[2] < M2.tuple[2] &&
		    M1.tuple[3] < M2.tuple[3] &&
		    M1.tuple[4] < M2.tuple[4] &&
		    M1.tuple[5] < M2.tuple[5] &&
		    M1.tuple[6] < M2.tuple[6] &&
		    M1.tuple[7] < M2.tuple[7] &&
		    M1.tuple[8] < M2.tuple[8] &&
		    M1.tuple[9] < M2.tuple[9] &&
		    M1.tuple[10] < M2.tuple[10] &&
		    M1.tuple[11] < M2.tuple[11] &&
		    M1.tuple[12] < M2.tuple[12] &&
		    M1.tuple[13] < M2.tuple[13] &&
		    M1.tuple[14] < M2.tuple[14] &&
		    M1.tuple[15] < M2.tuple[15]);
	}

	public static bool operator>(Matrix4d M1, Matrix4d M2) {
	    return (M1.tuple[0] > M2.tuple[0] &&
		    M1.tuple[1] > M2.tuple[1] &&
		    M1.tuple[2] > M2.tuple[2] &&
		    M1.tuple[3] > M2.tuple[3] &&
		    M1.tuple[4] > M2.tuple[4] &&
		    M1.tuple[5] > M2.tuple[5] &&
		    M1.tuple[6] > M2.tuple[6] &&
		    M1.tuple[7] > M2.tuple[7] &&
		    M1.tuple[8] > M2.tuple[8] &&
		    M1.tuple[9] > M2.tuple[9] &&
		    M1.tuple[10] > M2.tuple[10] &&
		    M1.tuple[11] > M2.tuple[11] &&
		    M1.tuple[12] > M2.tuple[12] &&
		    M1.tuple[13] > M2.tuple[13] &&
		    M1.tuple[14] > M2.tuple[14] &&
		    M1.tuple[15] > M2.tuple[15]);
	}

	public static bool operator<=(Matrix4d M1, Matrix4d M2) {
	    return (M1.tuple[0] <= M2.tuple[0] &&
		    M1.tuple[1] <= M2.tuple[1] &&
		    M1.tuple[2] <= M2.tuple[2] &&
		    M1.tuple[3] <= M2.tuple[3] &&
		    M1.tuple[4] <= M2.tuple[4] &&
		    M1.tuple[5] <= M2.tuple[5] &&
		    M1.tuple[6] <= M2.tuple[6] &&
		    M1.tuple[7] <= M2.tuple[7] &&
		    M1.tuple[8] <= M2.tuple[8] &&
		    M1.tuple[9] <= M2.tuple[9] &&
		    M1.tuple[10] <= M2.tuple[10] &&
		    M1.tuple[11] <= M2.tuple[11] &&
		    M1.tuple[12] <= M2.tuple[12] &&
		    M1.tuple[13] <= M2.tuple[13] &&
		    M1.tuple[14] <= M2.tuple[14] &&
		    M1.tuple[15] <= M2.tuple[15]);
	}

	public static bool operator>=(Matrix4d M1, Matrix4d M2) {
	    return (M1.tuple[0] >= M2.tuple[0] &&
		    M1.tuple[1] >= M2.tuple[1] &&
		    M1.tuple[2] >= M2.tuple[2] &&
		    M1.tuple[3] >= M2.tuple[3] &&
		    M1.tuple[4] >= M2.tuple[4] &&
		    M1.tuple[5] >= M2.tuple[5] &&
		    M1.tuple[6] >= M2.tuple[6] &&
		    M1.tuple[7] >= M2.tuple[7] &&
		    M1.tuple[8] >= M2.tuple[8] &&
		    M1.tuple[9] >= M2.tuple[9] &&
		    M1.tuple[10] >= M2.tuple[10] &&
		    M1.tuple[11] >= M2.tuple[11] &&
		    M1.tuple[12] >= M2.tuple[12] &&
		    M1.tuple[13] >= M2.tuple[13] &&
		    M1.tuple[14] >= M2.tuple[14] &&
		    M1.tuple[15] >= M2.tuple[15]);
	}

	public override bool Equals(object value) {
	    if(value == null) return false;
	    if (GetType() == value.GetType()) return this == (Matrix4d)value;
	    
	    return false;
	}

	public override int GetHashCode() {
	    return tuple.GetHashCode();
	}

	// operators
	public static Matrix4d operator+(Matrix4d M1, Matrix4d M2) {
	    return new Matrix4d(M1.tuple[0] + M2.tuple[0],
				M1.tuple[1] + M2.tuple[1],
				M1.tuple[2] + M2.tuple[2],
				M1.tuple[3] + M2.tuple[3],
				M1.tuple[4] + M2.tuple[4],
				M1.tuple[5] + M2.tuple[5],
				M1.tuple[6] + M2.tuple[6],
				M1.tuple[7] + M2.tuple[7],
				M1.tuple[8] + M2.tuple[8],
				M1.tuple[9] + M2.tuple[9],
				M1.tuple[10] + M2.tuple[10],
				M1.tuple[11] + M2.tuple[11],
				M1.tuple[12] + M2.tuple[12],
				M1.tuple[13] + M2.tuple[13],
				M1.tuple[14] + M2.tuple[14],
				M1.tuple[15] + M2.tuple[15]);
	}

	public static Matrix4d operator-(Matrix4d M1, Matrix4d M2) {
	    return new Matrix4d(M1.tuple[0] - M2.tuple[0],
				M1.tuple[1] - M2.tuple[1],
				M1.tuple[2] - M2.tuple[2],
				M1.tuple[3] - M2.tuple[3],
				M1.tuple[4] - M2.tuple[4],
				M1.tuple[5] - M2.tuple[5],
				M1.tuple[6] - M2.tuple[6],
				M1.tuple[7] - M2.tuple[7],
				M1.tuple[8] - M2.tuple[8],
				M1.tuple[9] - M2.tuple[9],
				M1.tuple[10] - M2.tuple[10],
				M1.tuple[11] - M2.tuple[11],
				M1.tuple[12] - M2.tuple[12],
				M1.tuple[13] - M2.tuple[13],
				M1.tuple[14] - M2.tuple[14],
				M1.tuple[15] - M2.tuple[15]);
	}

	public static Matrix4d operator*(Matrix4d M1, Matrix4d M2) {
	    return new Matrix4d(M1.tuple[0] * M2.tuple[0] + 
				M1.tuple[1] * M2.tuple[4] + 
				M1.tuple[2] * M2.tuple[8] +
				M1.tuple[3] * M2.tuple[12],

				M1.tuple[0] * M2.tuple[1] + 
				M1.tuple[1] * M2.tuple[5] + 
				M1.tuple[2] * M2.tuple[9] +
				M1.tuple[3] * M2.tuple[13],

				M1.tuple[0] * M2.tuple[2] + 
				M1.tuple[1] * M2.tuple[6] + 
				M1.tuple[2] * M2.tuple[10] +
				M1.tuple[3] * M2.tuple[14],

				M1.tuple[0] * M2.tuple[3] + 
				M1.tuple[1] * M2.tuple[7] + 
				M1.tuple[2] * M2.tuple[11] +
				M1.tuple[3] * M2.tuple[15],

				M1.tuple[4] * M2.tuple[0] + 
				M1.tuple[5] * M2.tuple[4] + 
				M1.tuple[6] * M2.tuple[8] +
				M1.tuple[7] * M2.tuple[12],

				M1.tuple[4] * M2.tuple[1] + 
				M1.tuple[5] * M2.tuple[5] + 
				M1.tuple[6] * M2.tuple[9] +
				M1.tuple[7] * M2.tuple[13],

				M1.tuple[4] * M2.tuple[2] + 
				M1.tuple[5] * M2.tuple[6] + 
				M1.tuple[6] * M2.tuple[10] +
				M1.tuple[7] * M2.tuple[14],

				M1.tuple[4] * M2.tuple[3] + 
				M1.tuple[5] * M2.tuple[7] + 
				M1.tuple[6] * M2.tuple[11] +
				M1.tuple[7] * M2.tuple[15],

				M1.tuple[8] * M2.tuple[0] + 
				M1.tuple[9] * M2.tuple[4] + 
				M1.tuple[10] * M2.tuple[8] +
				M1.tuple[11] * M2.tuple[12],

				M1.tuple[8] * M2.tuple[1] + 
				M1.tuple[9] * M2.tuple[5] + 
				M1.tuple[10] * M2.tuple[9] +
				M1.tuple[11] * M2.tuple[13],

				M1.tuple[8] * M2.tuple[2] + 
				M1.tuple[9] * M2.tuple[6] + 
				M1.tuple[10] * M2.tuple[10] +
				M1.tuple[11] * M2.tuple[14],

				M1.tuple[8] * M2.tuple[3] + 
				M1.tuple[9] * M2.tuple[7] + 
				M1.tuple[10] * M2.tuple[11] +
				M1.tuple[11] * M2.tuple[15],

				M1.tuple[12] * M2.tuple[0] + 
				M1.tuple[13] * M2.tuple[4] + 
				M1.tuple[14] * M2.tuple[8] +
				M1.tuple[15] * M2.tuple[12],

				M1.tuple[12] * M2.tuple[1] + 
				M1.tuple[13] * M2.tuple[5] + 
				M1.tuple[14] * M2.tuple[9] +
				M1.tuple[15] * M2.tuple[13],

				M1.tuple[12] * M2.tuple[2] + 
				M1.tuple[13] * M2.tuple[6] + 
				M1.tuple[14] * M2.tuple[10] +
				M1.tuple[15] * M2.tuple[14],

				M1.tuple[12] * M2.tuple[3] + 
				M1.tuple[13] * M2.tuple[7] + 
				M1.tuple[14] * M2.tuple[11] +
				M1.tuple[15] * M2.tuple[15]);
	}

	public static Matrix4d operator*(Matrix4d M, double Scalar) {
	    return new Matrix4d(M.tuple[0] * Scalar,
				M.tuple[1] * Scalar,
				M.tuple[2] * Scalar,
				M.tuple[3] * Scalar,
				M.tuple[4] * Scalar,
				M.tuple[5] * Scalar,
				M.tuple[6] * Scalar,
				M.tuple[7] * Scalar,
				M.tuple[8] * Scalar,
				M.tuple[9] * Scalar,
				M.tuple[10] * Scalar,
				M.tuple[11] * Scalar,
				M.tuple[12] * Scalar,
				M.tuple[13] * Scalar,
				M.tuple[14] * Scalar,
				M.tuple[15] * Scalar);
	}

	public static Matrix4d operator/(Matrix4d M, double Scalar) {
	    if(Scalar != 0.0f) {
		double reciproke = 1.0/Scalar;

		return new Matrix4d(M.tuple[0] * reciproke,
				    M.tuple[1] * reciproke,
				    M.tuple[2] * reciproke,
				    M.tuple[3] * reciproke,
				    M.tuple[4] * reciproke,
				    M.tuple[5] * reciproke,
				    M.tuple[6] * reciproke,
				    M.tuple[7] * reciproke,
				    M.tuple[8] * reciproke,
				    M.tuple[9] * reciproke,
				    M.tuple[10] * reciproke,
				    M.tuple[11] * reciproke,
				    M.tuple[12] * reciproke,
				    M.tuple[13] * reciproke,
				    M.tuple[14] * reciproke,
				    M.tuple[15] * reciproke);
	    }
	    else {
		return new Matrix4d(double.MaxValue,
				    double.MaxValue,
				    double.MaxValue,
				    double.MaxValue,
				    double.MaxValue,
				    double.MaxValue,
				    double.MaxValue,
				    double.MaxValue,
				    double.MaxValue,
				    double.MaxValue,
				    double.MaxValue,
				    double.MaxValue,
				    double.MaxValue,
				    double.MaxValue,
				    double.MaxValue,
				    double.MaxValue);
	    }
	}

	public static Vector4d operator*(Matrix4d M, Vector4d V) {
	    return new Vector4d(M.tuple[0] * V.X + 
				M.tuple[1] * V.Y + 
				M.tuple[2] * V.Z +
				M.tuple[3] * V.W,
				
				M.tuple[4] * V.X + 
				M.tuple[5] * V.Y + 
				M.tuple[6] * V.Z +
				M.tuple[7] * V.W,
				
				M.tuple[8] * V.X + 
				M.tuple[9] * V.Y + 
				M.tuple[10] * V.Z +
				M.tuple[11] * V.W,

				M.tuple[12] * V.X + 
				M.tuple[13] * V.Y + 
				M.tuple[14] * V.Z +
				M.tuple[15] * V.W);
	}

	public static Matrix4d operator-(Matrix4d M) {
	    return new Matrix4d(-M.tuple[0],
				-M.tuple[1],
				-M.tuple[2],
				-M.tuple[3],
				-M.tuple[4],
				-M.tuple[5],
				-M.tuple[6],
				-M.tuple[7],
				-M.tuple[8],
				-M.tuple[9],
				-M.tuple[10],
				-M.tuple[11],
				-M.tuple[12],
				-M.tuple[13],
				-M.tuple[14],
				-M.tuple[15]);
	}

	public Matrix4d Transpose() {
	    return new Matrix4d(tuple[0],
				tuple[4],
				tuple[8],
				tuple[12],
				tuple[1],
				tuple[5],
				tuple[9],
				tuple[13],
				tuple[2],
				tuple[6],
				tuple[10],
				tuple[14],
				tuple[3],
				tuple[7],
				tuple[11],
				tuple[15]);
	}

	private double[] tuple;
    }

}
